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Abstract. Shrock and Wu have given numerical values for the exponential 
' growth rate of the number of spanning trees in Euclidean lattices. We give a 

new technique for numerical evaluation that gives much more precise values, 
together with rigorous bounds on the accuracy. In particular, the new values 
\Q ' resolve one of their questions. 
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§1. Introduction. 

o 
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Since the time of Kirchhoff ( |1847|) , physicists and mathematicians have been inter- 
ested in enumerating spanning trees. One aspect of this endeavor has been to evaluate 
or to estimate asymptotics of the growth rate of the number of spanning trees in large 
graphs. Additional interest in the asymptotic growth rate arises in ergodic theory, since 



the exponential rate is also the entropy of a natural and important system, the so-called 
j> | uniform spanning forest. See Pemantle (|1991| ), Burton and Pemantle ( |1993|) , Benjamini, 



Lyons, Peres, and Schramm (|2001|) , Lyons ( |1998| ), and Lyons (|2003| ) for explanations and 
information about the uniform spanning forest. Some modern enumeration efforts include 



Burton and Pemantle (|1993| ), Shrock and Wu (|2000|) , Lyons ( |2003| ); see also the references 



therein. In many cases, one can express the main term of the asymptotics by an integral 
formula. For example, if t(G) denotes the number of spanning trees of a graph G and if 
G n are the graphs induced by cubes of side length n in the hypercubic lattice Z d , then 
it is well known (and rederived in both Burton and Pemantle ( |1993| ) and Shrock and Wu 
fl2000|) ) that the thermodynamic limit 

hd ■= lim logT(G n ) 

n— »oo fl 
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can be expressed as 



hd= f log ( 2d — 2 cos 2nXi ] dx = log(2d)+[ log [ 1 — — cos 27TXj | da;. (1.1) 

It is also well known (due to its connection with the dimer problem) that hi = 4G/7T, where 
G := Efclo( -1 ) fc /( 2/c + l f is Catalan's constant (see, e.g., Kasteleyn (|TMT| ) or Montroll 
( p.964[ )). No values of ^ for any d > 3 are known in simple terms of other known constants 
and functions. Shrock and Wu ( |2UUU| ) evaluated these integrals in higher dimensions by 
numerical methods and found one particularly intriguing value: /14 = 2.0000(5). They 
suggested that /14 may be exactly 2, which would be quite surprising. Indeed, it would 
be extraordinary for a natural system without parameters to have a natural-log entropy 
that is a non-zero integer. This would have been the first such example to our knowledge. 
However, we shall see that h& is extraordinarily close to 2, but not, in fact, exactly 2. We 
shall also give more accurate values for other hd with rigorous bounds on their accuracy. 

The numerical evaluation of hd is problematic if one wants to use the formula |( 1 . 1 )| , due 
to the difficulty of accurate integration in higher dimensions. Therefore, Shrock and Wu 
( |2000p gave an interesting large-d asymptotic expansion of hd to order 1/d 6 ; we have given 



more terms below, to show that not all coefficients are positive, as one might otherwise 
believe, and to illustrate that this is indeed only an asymptotic expansion, not a convergent 
series: 



h d = log(2d) 



1 3 7 45 269 805 3615 

+ ~T7> To + 7T7, — 77 + 7777^7 + TTTTT^F + 7777^ + 



Ad 16 d 2 32 d 3 128 d 4 384 d 5 512 d 6 1024 d 7 
23205 144963 2187031 40225409 1277353077 



+ 



4096 d 8 10240 d 9 8192 d 10 16384 c/ 11 65536 d 12 
66817216455 271891453119 / 1 

+ o 



458752 d 13 262144 d 14 \d 15 

However, it is difficult to know how many terms of this expansion to use; Shrock and Wu 
( |2000|) used this series to report h$ = 2.243 and h§ = 2.437. For smaller d, Shrock and Wu 



( [2000D used numerical integration to find that h 3 = 1.6741481(1) and that h 4 = 2.0000(5). 
However, the accuracy of h 3 is off by several orders of magnitude. 

To obtain greater accuracy and enable us to prove that h± ^ 2, we shall use a new 
formula, namely, 

oo 

h d = log(2d)-Y,Pd(k)/k, (1.2) 

fc=i 
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where Pd(k) is the probability that simple nearest-neighbor random walk on the hypercubic 
lattice Z d returns to its starting point after k steps. (A generalization of this formula is 
due to Lyons (|2003| ) . ) Even this formula is slightly problematic to use, since Pd(k) is the 
sum of a large number of binomial coefficients for large k. Because of the large number of 
small terms, it is important to calculate the sum as an exact rational before converting to 
a real approximation. Fortunately, there is a simple recursion formula that enables quicker 
computation. In addition, we explain how to estimate the tail of the series in |(1.2j. 



In the remainder of this note, we first state our numerical results and then derive 
the simple but crucial |(1.2)| . Next, we explain how to compute Pd(k) quickly and how to 
approximate the error, and finally prove rigorous bounds. We shall also briefly discuss 
body-centered cubic lattices. We end by discussing an alternative approach that was 
brought to our attention after a first version of this article was submitted. 



§2. Results. 

The numerical results are, to an accuracy we believe includes all reported digits, in 
the following table: 



d 


h d 




d 




3 


1.67338930297 




12 


3.1557714292824 


4 


1.999707644517 




13 


3.2376421551842 


5 


2.2424880598113 




14 


3.31330031802725 


6 


2.43662696200071 




15 


3.383624540390254 


7 


2.5986763042 




16 


3.449318935201 


8 


2.73786766385 




17 


3.510956551787645 


9 


2.859910142340 




18 


3.56901006528479 


10 


2.968594484443 




19 


3.62387396384455 


11 


3.066571824248 




20 


3.67588091671257 



These arise as follows. To prove |(1.2)| , use the Maclaurin series for log(l — z) to find 

f, / \ ^ i r h d l fc 

/ log 1 > cos Itxx j \ dx = — y — I — > cos27rx,- ax. 

= -E 



k=i 
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Clearly, we have Pd(k) = for k odd and p\{2k) = ( 2 fc fc )/2 2fc . It is well known that 
P2(2k) = pi(2k) 2 (e.g., one step of a random walk in Z 2 can be made by taking one step 
in each of the directions ±(l/\/2, l/\/2) and ±(l/v2, — 1/\/2), independently). For fast 
computation of other return probabilities, write /(<i, k) := [2d) 2k pd(2k) for the number of 
nearest- neighbor walks of length 2k in Z d that start and end at the origin. Such a walk 
has the property that its projection to the first d\ coordinates also starts and ends at 
the origin, while the number of steps in the first d\ directions may be any even number 
between and 2k. Similar reasoning shows that 



k 

f(di + d 2 , k) = J2 (l ) f(di, r)f(d 2 , k - r) 

r=0 ^ ' 



This allows one to reduce the computation of Pd(') to the values of P[d/2] (•) an d P\d/2~\ (•)• 
Once we have these values, it is simple to compute partial sums for |(1.2) . Since all 



terms of the series are positive, each such partial sum gives a rigorous upper bound for the 



true value of hd- Merely summing the first 13 terms of the series in |(1.2)| for d = 4 gives a 
rigorous proof that /14 < 2. To get an lower bound for hd, it suffices to bound above the 
remainder. It is well known (see, e.g., Spitzer ( |1976|) , Section 7) that 



d 



d/2 



A more accurate approximation is 

as shown by Ball and Sterbenz (|2003| ). As this suggests, we believe that the right-hand 
side of |(2.1)| is actually larger than the left-hand side; indeed, this appears to be true for 



all d and k, not merely for large k, though it has been proved only for large k and small 
d. That is, we have 

d \ d/2 
■J 



M2k) ~ 2 <4^ J (2 ' 3) 



for all k when 1 < d < 6 and for all large k (if not all k) when d > 7; see Ball and Sterbenz 
( |2003| ). Since the sum over k > r, any r > 0, of the right-hand sides of either |(2.1)| or |(2.2 



can be expressed via the Hurwitz zeta function, for which Euler-Maclaurin summation 
approximations are readily available, we obtain the very accurate values reported in the 
tables above by summing relatively few terms. Excellent accuracy is already available after 
just 10 terms, but we have used 1000 terms for 3 < d < 6, 100 terms for 7 < d < 10, and 
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80 terms for 11 < d < 20. In addition, by using [2.3) and partial sums of 1000 terms, we 
get the rigorous bounds 

1.6733893024176978 <h 3 < 1.6733917596720884 
1.9997076445004571 <h 4 < 1.9997076951104138 
2.242488059810819 <h 5 < 2.2424880610724065 
2.436626962000695 <h 6 < 2.4366269620369234 . 

One can use similar estimates to improve the accuracy of the asymptotics for body- 
centered hypercubic lattices. As shown by Shrock and Wu ( |2000| ), the exponential growth 
rate of the number of spanning trees in d dimensions is 

1 00 1 

h h d cc = dlog2--J2 l Pi(k) d . (2.4) 

k=i 



It is straightforward to show that 



1 ! 



by use of Stirling's approximation. We sum 1000 terms of |(2.4)| and bound the remainder. 



Using |(2.3)| for a lower bound and |(2.5)| for an upper bound, we find 



1.9901914178466 <h% cc < 1.9901914178472 
2.732957535468933 <h\ cc < 2.7329575354689455 



These agree with the estimates of Shrock and Wu ( |2000p , but give about 3 times as many 
digits. 

After a first version of this article was submitted for publication, Alan Sokal kindly 
brought to our attention some related calculations by Sokal and Starinets ( |2001|) . Up to 
a constant, the entropy studied here is equal to the free energy g^il/d) studied there 
(see equation (A. 2) of their paper). Their formula (A. 6) shows, then, that 

f°° e~ l 

h d = log(2d) + / — [1 - h{t/df] dt , (2.6) 
Jo t 

where Iq is the modified Bessel function. In this way, hd can be estimated by numerical 
integration in only one dimension, which can be accomplished very quickly. The disad- 
vantage, however, is that the integrand decays rather slowly. As noted in Appendix A. 2 
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of Sokal and Starinets ( [2001 ), one can improve the precision dramatically by numerical 
integration up to some cut-off, then symbolic integration of the tail with an asymptotic 
formula replacing Iq. Even so, not all the numerical values reported in Sokal and Starinets 
( POQI ) are correct in all their digits, as can be seen by comparison with our tables and 
our rigorous bounds. The second disadvantage of |(2.6)| is that it is less straightforward 
to provide rigorous bounds. For this purpose, one has to treat carefully the technique of 
numerical integration, as well as evaluation and bounding of Iq. Some of this is discussed 
in Appendix B of Hara and Slade ( |1992j ). By comparison, our technique requires only the 
calculation of rational numbers, as well as one simple logarithm calculation. 
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